Time-Dependent and Steady- State Gutz wilier approach for nonequilibrium 

transport in nanostructures 
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We extend the time-dependent Gutzwiller variational approach, recently introduced by Schiro 
and Fabrizio, Phys. Rev. Lett. 105 076401 (2010), to impurity problems. Furthermore, we derive 
a consistent theory for the steady state, and show its equivalence with the previously introduced 
nonequilibrium steady-state extension of the Gutzwiller approach. The method is shown to be able 
to capture dissipation in the leads, so that a steady state is reached after a sufficiently long relaxation 
time. The time-dependent method is applied to the single orbital Anderson impurity model at half- 
filling, modeling a quantum dot coupled to two leads. In these first exploratory calculations the 
Gutzwiller projector is limited to act only on the impurity. The strengths and the limitations of this 
approximation are assessed via comparison with state of the art continuous time quantum Monte 
Carlo results. Finally, we discuss how the method can be systematically improved by extending the 
region of action of the Gutzwiller projector. 

PACS numbers: 73.63.Kv, 73.63.-b, 71.27.+a 



I. INTRODUCTION 

The impressive advances in nanoelectronics has en- 
abled contacting of very small objects, such as quantum 
dots, molecular junctions and nanowires, with metal- 
lic electrodes, enabling very accurate measuraments of 
the source-drain current counting individual electrons as 
they tunnel across the contact»i~— Because of the low di- 
mensionality of the contact region correlations grow in 
strength and may stabilize a local magnetic moment that 
influences electron tunneling^ 

Due to the interplay between strong electron cor- 
relation and out-of-cquilibrium effects, the theoret- 
ical study of these systems is extremely compli- 
cated. Apart from many-body Keldysh pertur- 
bation theory— >S many innovative approaches have 
been explored^ such as diagrammatic quantum Monte 
Carlo on the Keldysh contour^— field theory tech- 
niques^^ time dependent density matrix renormaliza- 
tion group ; 14 ' 15 flow equation method ; 16 : 17 functional 
renormalization group ; 18 ' 19 perturbative renormalization 
groups 0- — master equations^! iterative path integral 
approaches^ strong-coupling expansions r§r— real time 
numerical renormalization group ; 30 : 31 scattering Bethe 
Ansatzj 3 - 2 - and imaginary-time nonequilibrium quantum 
Monte Carlo^i 

Unfortunately, most of these techniques arc computa- 
tionally very demanding, limiting their application only 
to simple models. For realistic cases, e.g., tunneling 
across a molecule or a transition metal atom — where 
many orbitals participate to magnetism and affect con- 
ductance — these methods become generally intractable. 
In this general context a sufficiently simple approach to 
deal with realistic systems, even one being less accurate 
than those previously mentioned, would be extremely 
useful. 

With this motivation, an extension of the Gutzwiller 



variational metho d ' to nonequilibrium steady-state 
transport was recently proposed in Ref. H3, and applied 
to the single-orbital Anderson impurity model at half- 
filling. In parallel, a generalization of the time- dependent 
Gutzwiller approach 3 ^ was developed by Schiro and Fab- 
rizio^ 9 - on the basis of the Dirac-Frenkcl variational prin- 
ciple^2r— and applied to study the quantum dynamics 
induced by an interaction quench in the single band Hub- 
bard model. 

The first achievement of this paper is to generalize 
the time dependent Gutzwiller approach of Ref. |H to 
impurity problems and derive a consistent theory for 
the steady state. We find that dissipation occurs en- 
tirely within the leads as expected^ and the steady-state 
reached after a sufficiently long thermalization time is the 
same as the one previously derived in Ref. H3- We believe 
that our scheme represents a very interesting step toward 
the description of steady state non-equilibrium problems. 
In fact, with this genuine variational principle one can, 
in principle, always obtain more accurate results by ex- 
tending the variational freedom of the wavcf unction. 

The second achievement of this paper is to investigate 
the reliability of the Gutzwiller variational ansatz in the 
case of a Gutzwiller projector acting only on the impurity 
(as was assumed in Ref. H3) , which is the simplest con- 
ceivable Gutzwiller variational function for an impurity 
model. In this investigation we employ the single-orbital 
Anderson impurity model as a prototype, and study two 
different ways to prepare the initial state before start- 
ing the time evolution: the "interaction quench" and the 
"bias quench" . We compare our calculations with the 
quantum Monte Carlo results of Ref. [Hi The compari- 
son shows that our results are quantitatively predictive at 
short times. At longer times, instead, the Gutzwiller dy- 
namics can develop unphysical features, due to the over- 
simplified variational wavefunction. We expect that the 
quality of the result would be considerably improved if 
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the Gutzwillcr projector acted also on a portion of leads. 
In particular, the resulting additional variational freedom 
would allow us to better account for the Kondo effect. 

The outline of this paper is as follows. In Sec. [TT] the 
time dependent Gutzwillcr method is generalized to im- 
purity models on the basis of the Dirac-Frenkel varia- 
tional principle, and specialized to the single orbital An- 
derson impurity model at half-filling. In Sec. lIIII the vari- 
ational approach for the steady state is derived on the 
same basis. Furthermore, several important conceptual 
differences with respect to the theory at equilibrium are 
discussed. In Sec. IIVI the Gutzwillcr variational results 
for the interaction and the bias quench are presented. 
Both the transient and the long time dynamics are ana- 
lyzed. Finally, Sec. [V] is devoted to conclusions. 

II. TIME DEPENDENT GUTZWILLER 
METHOD FOR IMPURITY MODELS 

The Dirac-Frenkel variational principle ; 40 ' 41 identifies 
the Schrodinger quantum dynamics of the system 

|*(t)> = e-^i) (1) 
with the saddle point of the action 

5[*(f)] = f dr£[tt(r)] (2) 
£[*(r)] = (*(r)|*a r -'W|*(r)>, (3) 

searched in the set of all the possible time evolutions 
|\D'(<)) in the Hilbcrt space — with fixed boundaries at tj 
and tf. The advantage of the Dirac-Frenkel formulation 
of the dynamics is that it allows us to build up a well 
founded variational approximation scheme. In fact, the 
saddle point of S can be searched in a proper subclass of 
the Hilbert space, chosen on the basis of physical moti- 
vations related with the specific system of interest. An 
important example is the time- dependent Hartree-Fock 
method, that can be derived as the saddle point of the 
Dirac-Frenkel action S assuming that |^) spans the set 
of all Slater determinants. 

As shown in Ref. |42| by Schiro and Fabrizio, the Dirac- 
Frcnkcl principle provides a solid basis also to the time- 
dependent Gutzwillcr variational method, that was pre- 
viously introduced by the same authors in Ref. and 
used to study the quantum dynamics induced by an in- 
teraction quench in the single band Hubbard model. The 
time-dependent Gutzwillcr dynamics is expected to im- 
prove considerably the time-dependent Hartree-Fock dy- 
namics for interacting systems, as the Gutzwillcr varia- 
tional space enlarges the set of Slater determinants by 
means of the Gutzwiller projector, that is able to modify 
the weight of local electronic configurations in accordance 
with the interaction. In this section we show how to gen- 
eralize the time-dependent Gutzwiller scheme to study 
impurity problems. 



A. Gutzwiller dynamics of the single-orbital 
Anderson model 

For simplicity we employ the single-orbital Anderson 
model at half filling — modeling a quantum dot coupled 
to two leads — as a prototype for this scheme. The 
corresponding Hamiltonian is 

n = t + V+U, (4) 

consisting of the following terms 

aka 

^ = E^ C ^X + H.c. (6) 

aka 

U = \in A - l) 2 , (7) 

where , creates a conduction electron on the left a = 
— 1 or right a = 1 lead with quantum number k and spin 
a while d\, creates an electron in the dot with spin <j, and 
fi is the quantization volume of the system. 

If the initial state |^,) is particle- hole symmetric then 
so is also ^(t)) at any time. This allows us to restrict 
our attention to the particlc-holc symmetric variational 
space. We make the subsequent variational ansatz for 
the time-dependent wavefunction^ 

|*> = Vg\*o) = e- iS «V d \* ) , (8) 
S d = Q2 P Q2 (9) 

y d = ^02^02 + y / 2 _ ^A (io) 

where Z02 and #02 are real parameters, 

P 02 = | t 4)(t4|+|o)(o|, P, = |tXt|+|4>U| (11) 

are impurity projection operators, and \^o) is a Slater 
determinant. It can be easily verified that P02 = 
(#|P 2 I*) = Iq 2 /2, so that P 02 and 6» 02 can be used 
as variational parameters. Furthermore, it can be easily 
shown that the wavefunction defined above satisfies 
the Gutzwillcr conditions 

(* |^P G |*o) = 1 (12) 
(*o|P^ G 4^l*o) = \ (13) 

for all possible values of < P02 < 1 and #02 ■ Notice 
that P02 is just twice the expectation value of the double- 
occupancy of the dot. 

Using Eqs. (fl"2"]) and (fT5|) one finds that the average 
energy £ with respect to the Gutzwiller variational wave- 
function [Eq. ©] is given by 

£[P 2,#02,*0] = (*|#|«) 

= (*o|^|*o) + C/Po2/2, (14) 
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where 



H$ = t + RV 



R = 2^02(1-^02)008(002). 



(15) 
(16) 



Furthermore, from Eqs. (|T2"j) and (|T3"|) we obtain the fol- 
lowing analytical expression for the Lagrangian 4 ^ 



In the particle-hole symmetric case the left and right cur- 
rents are equal in magnitude and the lead index a can 
be dropped, I = |7 a |. 



III. GUTZWILLER THEORY FOR THE 
STEADY STATE 



02 j o 02 j ^0 



- £[P02, 002,^0}- 



The corresponding equations of motion are 
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(17) 

(18) 
(19) 
(20) 



Let us now rewrite Eqs. (|18 p -(|2U p in a soluble form. 
Substituting Eqs. (0 and (15]) in Eqs. §Tg\) and ([19)) 

gives 



BR 



112 



at 
ag 02 
at 



-V(t) 



_9_R 



[/ , BR 
2 +V ^BW 2 



(21) 
(22) 



where V(t) = (* (*)| V |#o(*))- The time evolution of 
I "Jo), see Eq. pO|. can be formulated in the Heisenberg 
picture as follows: 



at 

RV{dUa) - 



dt 



'aka a I 



^Rvn{cl ka c a , k , a ) 



a' k' 



(23) 



|(e fc - efe')°(4/a'fev)+ 



W (Vn^c^J - ^<4X>) } , (24) 

where, for simplicity, a constant Vfe = V has been 
assumed, and the averages are taken with respect to 
|^o(*))- Notice that the occupancy of the dot {d\d a ) 
is 1/2 due to particle-hole symmetry. 

The system defined by Eqs. (|2ip - (|2"4"|) constitutes a set 
of coupled first order differential equations, which can be 
readily solved numerically. From the time evolution of 
the parameters that appear in Eqs. ([2"T j) - (f2"4")) is possible 
to calculate the expectation value of any desired observ- 
able. In particular, the expectation value of the current 
is obtained as 



■ V ST 



ka 



ka 



ko&a) 



(25) 



Let us consider a quantum system whose dynamics is 
generated by the action [Eq. ((2J)] in a Hilbert space X. 
The condition of stationarity for a state \^>) is that 



6^ 



V|<5*)i_|*> ex. 



(26) 



a t *=o 



In this section we will derive the Gutzwillcr theory for the 
steady state starting from the general condition [Eq. (|2"6"j) ] 
applied to the Dirac-Frenkel dynamics in the Gutzwiller 
variational space. 

For clarity we will focus on the single orbital Anderson 
impurity model at half- filling and infinite leads. This will 
facilitate the subsequent discussion of the most general 
case. Furthermore, we will assume — without any loss of 
generality 4 ^ — a one-dimensional representation of the 
Hamiltonian of the leads [Eq. ([5])] 



"T - tRR ' C ^Ra C a R'a > 



(27) 



aa RR' 



where the site labels R and R' run over all integer num- 
bers from 1 to 00. 



A. Infinite systems and inequivalent 
representations 

Before proceeding with the derivation of the Gutzwiller 
scheme, we need a preliminary discussion of the problem 
of the steady state in relation with the general structure 
of infinite fermionic systems. 

Let us define the C* -algebr a 45 ' 46 A generated by the 
local fermionic ladder operators c Rcr — satisfying the 
canonical anticommutation relations. 

We associate to a generic state \^f) of the Anderson 
system [Eq. Q] its corresponding sector X\$>\, defined as 
the Hilbert space spanned by all the states ^4| V I / ), with 
A e A. The linear space X[^>] is the basis of a represen- 
tation of A, i.e., 



G A\<t>) e J[tf] VAeA. 



(28) 



The physical meaning of Eq. (f2"8")) is that X[*5f] repre- 
sents a class of "macroscopically equivalent" states, that 
differs from \*5f) only by local modifications. In this gen- 
eralized sense, a sector can be regarded as a phase of the 
system^ 

Note that, in general, it is not guaranteed that the 
sectors of two different states coincide. In fact, the Von 
Neumann theorem docs not apply to infinite systems^ 5 - 
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and the obscrvablcs generally admit a big variety of in- 
equivalent representations. 

Let us define the initial state |W; n ) of two unconnected 
leads prepared at chemical potentials \x a = a$/2, where 
$ corresponds to the applied bias. Once the tunneling 
region [Eq. ((5])] is connected to the leads, the initial state 
| "fin) evolves according to the Schrodinger equation 



id t \*(t)) = H\*(t)) , 



(29) 



with £ I[^in] for any finite time t. A very im- 

portant point for our subsequent discussion is that the 
steady state l^s), that is reached after an infinite time, 
does not belong to Z^iJ. 

Let us prove our statement. It is clear that all the 
vectors that belong to the same sector share the same 
average left-current, which is defined as follows: 



^[*]=lj?lTr<*lEE / *»l*>' 

a ReA 



(30) 



where the symbol A j" denotes a sequence of bounded 
subregions A of the system that increase to infinity, |A| 
is the number of sites in A, 



iRa 



i[U,N J 



Rati 



(31) 



and Nn a represents the total number of electrons to the 
left of the site identified by R and a. In fact, the limit in 
Eq. (f3"0")) is not affected by the contribution of any local 
operator A, i.e., 



J[tf ] = I[AV] vi g A ■ 



(32) 



The meaning of Eq. ([3"2")) is that two macroscopically 
equivalent states share the same average current. In this 
generalized sensed the average current can be regarded 
as an "order parameter", that enables us to distinguish 
different sectors (phases) of states. 

This simple observation makes it easy to prove the 
statement that \^s) does not belong to I[^j n ]. In fact, 
after an infinite time, all expectation values of the local 
current in Eq. (|30[) are non-zero, so that 



I[^S] ± I[*in] = . 



(33) 



The thesis follows by contradiction from the fact proven 
above that if l^m) and \^s) would belong to the same 
sector then they would share the same average current. 
In particular, this observation implies that the limit 



I* 



3) = lim 



(34) 



can not be interpreted as a limit in the norm induced by 
the scalar product in a Hilbert space, but is well defined 
only in the following "weak" sense 

(9 a \A\*s) = lim(*(i)|i|*(t)) WAeA. (35) 

t— >oo 

Furthermore, the same argument shows that steady 
states of Hamiltonians that differ only locally — e.g., 



at the impurity — belong, in general, to inequivalent 
representations of the obscrvables. In fact, the steady 
state average current depends on the local interaction at 
the impurity. 

Notice that the argument above does not apply to the 
equilibrium the average current is zero indepen- 

dently on the junction at zero bias. In fact, in this case 
(and only in this case), |\&i n ) is the ground state of T, 
so that the action to connect the tunneling region to the 
leads corresponds to perturb locally the ground state of 
T with 



SH = V+U. 



(36) 



see Eqs. Q-Q. This implies that |^ s > is the ground 
state of H, which belongs to the same sector of |\&i n )- 
This is ensured by well established general results con- 
cerning the "stability of thermal sectors under local per- 
turbations" , see Ref. l46l 

In conclusion, the steady states of biased systems with 
different junctions belong to disjoint representations of 
the observables. This is a fundamental difference between 
equilibrium and out-of-equilibrium problems, and will be 
central in our subsequent discussion of the Gutzwiller 
theory for the nonequilibrium steady state. 



B. Steady state theory for the single orbital 
Anderson impurity model at half-filling 

In this section we show how the Gutzwiller steady-state 
theory for the single orbital Anderson impurity model 
at half-fillingS can be derived by means of the general 
scheme introduced above. 

From Eq. (JTTJ) we have that the stationarity condition 
[Eq. ([2"o]) ] is equivalent to 



5£(^ ,P 02 ,9 02 ) 

d£(y ,P 02 ,6 02 ) 
dP 02 

d£(^o,p Q2 ,e 02 ) 
de 02 



(37) 
(38) 
(39) 



where £ is given by Eq. (fT4"]) and \5^o) G Z^o] is orthog- 
onal to |* ). 

Let us assume that the Gutzwiller parameters Pq 2 and 
#02 converge to proper steady state values after a relax- 
ation time T due to dissipation within the leads ^ see 



Sec. lIVI This implies that for t > T only the Slater deter- 
minant j^o ) evolves in time, and its dynamics is induced 
by the effective steady-state Lagrangian 



(40) 



£o[*o] = (*o\idt-K\* ) 



where R is obtained from Pq 2 and 9q 2 using Eq. (jT6 
Equivalently, for t > T 



|*oW> = e-^ (t - T) |*o(T)). 



(41) 
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From (i) the uniqueness of the steady state and (ii) the 
observation that |*o(T)) and |*o(0)) differ only locally 
we obtain that, after an infinite time, the Slater determi- 
nant is given by 



K) = Urn e^" 4 |v!/ (0)), 



(42) 



which is simply the steady state of Hq . The above ar- 
gument reproduces the variational ansatz for the steady 
state used in Ref. [13 from a very clear perspective. In 
order to not evolve in time with Eqs. (fT8|) - ([20l) . the Slater 
determinant of a variational state \^>) of the form 
[Eq. flSJ)] is necessarily the steady state \ ^q ) of Hq, where 
R is obtained from P 02 and 802 using Eq. (TIB)) . 

This consideration facilitates the search for the steady 
state. In fact, we can restrict the search of the state 
that satisfy the stationarity conditions [Eqs. (j37|) - ([39)) ] to 
the variational subset of "candidate" Gutzwiller steady 
states of the form defined above, 



m = v G m 



(43) 



each one belonging to a distinct representation of the ob- 
servables. In our case, this amounts to search for the 
solution of Eqs. ((37)) and ([38]) only, as Eq. (f39| is auto- 
matically satisfied by From the stationarity condi- 
tions [Eqs. (|3T|) and (|38p)] we obtain that 



(*?ivi*?; 



dR 




0. 



7(12 



(44) 
(45) 



which ensure that P02 and #02 do not evolve in time, 
see Eqs. (|2TJ and p2]l). Notice that V |^) can be 
readily calculated by means of the Keldysh formalism^ or 
scattering theory ; 47 ' 48 as Hq is a quadratic Hamiltonian, 
see section llV Al 

Eq. (|4"S")) implies that the Gutzwiller projector of the 
Gutzwiller steady state is real, as it implies that #02 is a 
multiple of it, see Eq. ((9]). Eq. (|44[) allows us to calculate 
the steady state value of P02, as R is a known function 
of P 02 , see Eq. (TTB]). 

In conclusion, we have shown that if a Gutzwiller 
steady-state exists then it is determined by Eqs. (|43p- 

(SSJ). 

Note that, from the general point of view of Sec. IHI Al 
Eqs. (|3"7} - ([3"9"| establish the condition of stationarity of 
a given variational state I*) in its own sector I[$>] = 
I^q ], which is indeed, contrarily to the equilibrium case, 
a state- dependent variational space. 



Steady state theory for multi-orbital Anderson 
Models 



In this section we generalize the Gutzwiller steady- 
state theory to a general impurity model representing two 



infinite leads connected through a generic multi-orbital 
junction 

H = f + V+U, (46) 
where 14 is the Hamiltonian of the junction, and 

f = Y.^ c lk,c aka (47) 



aka n 



Vn 



(48) 



where c aka creates a conduction electron on the left a = 
— 1 or right a = 1 lead with quantum number k and spin 
a, (fi na creates an electron in the junction with quantum 
label n and spin cr, and f2 is the quantization volume of 
the system. 

We employ the general </>-matrix formalism, see Ref. [4^ 
and references therein, assuming, for simplicity, a spin 
rotationally invariant Gutzwiller variational state. Let 
us introduce the so-called natural-basis^ operators f na , 
i.e., the ladder operators such that 

(*o\flf m < 7 \*o)=6 nm n n Vn,m (49) 

where < n° < 1 are the eigenvalues of the local density 
matrix 



(*o|4,CI*o) = p1 



(50) 



Furthermore, we introduce the matrix representation of 
the operators f na 



(Fj)., = (r.ljur 



j 1 ■ 



(51) 



where \Ti) are many-body Fock states expressed in the 
/ ncr -basis. 

It can be shown that the Gutzwiller variational func- 
tion is parametrized by cf> and l^o) satisfying the 
Gutzwiller constraints 



(52) 
(53) 



Tr ((/>+</>) = 1 , 
Tr (^4>F^F m(r ) =n°J nm , 

and the variational energy is given by 

f[0,*o] = (*o|#ol*o) +Tr(0t[/0) , (54) 

where 

R" = t + V R (55) 
Vn = EE ^%^ /^ fcCT (56) 

aka n * 

Tr((/) t F^ 6F ) 

(57) 



R 



Ui 



(58) 
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It can be easily shown that the Dirac-Frenkcl La- 
grangian is given by 

£[0,*o,A] - (tf |*aetto)-£[&tfo] 

+ (^o\V G (d t V G ) |* ) + C[X; cf>, * ] , (59) 

where C is a Lagrange multiplier term that ensures that 
the Gutzwiller constraints - which can be regarded as 
holonomic constraints - are satisfied at any time. 

Let us assume to have solved the corresponding La- 
grange equations. Similarly to the case of the single or- 
bital Anderson impurity model, the matrix <j), the varia- 
tional density matrix n° and the Lagrange multipliers A 
converge to proper steady state values after a transient 
time T. This implies that for t > T only the Slater de- 
terminant |\&o) evolves in time by means of the dynamics 
induced by the effective steady-state Lagrangian 

£ [*0] = (* | i0 t -[W O R +E fUma] 1*0) , (60) 

nra cr 

where R is obtained from <j> and n° using Eq. (|58[) . In 
order to not evolve in time, \^o) is necessarily the steady 
state of 

HI [R, A] EE n§ + £ A„ m £ fnafma , (61) 
nrn a 

where R is obtained from <p and n° using Eq. (|58p and A 
ensures that |^o) satisfies Eq. (|49|) . 

From Eq. we deduce that the stationarity condi- 
tion [Eq. ((26|) ] is equivalent to 

<5£(^ o ,0) = (62) 

for all the variations of <fi and \5^o) G ^[^o] such that 
l^^o) _L |* ) and the "holonomic" constraints [Eqs. ([52]) 
and (j5"3")) ] are satisfied. Similarly to the case of the single 
orbital Anderson impurity model, the argument above 
facilitates the search of the steady state, as the solution 
of Eq. (|62p can be searched in the restricted variational 
subspace for the nonequilibrium steady state. 

We point out that Eq. (|6"2")l establishes the condition 
of stationarity of a variational state in that is 

a state- dependent variational space. 

Notice that, in principle, the general scheme described 
above allows us to study any impurity system with an 
arbitrary level of accuracy. Let us consider, for instance, 
the single orbital Anderson impurity model. One can 
treat the impurity and a portion of leads as a multi- 
orbital junction and study the system with the general 
method described above. The resulting increased vari- 
ational freedom would take into account also the local 
magnetic correlations between the impurity and the sur- 
rounding lead electrons. The so obtained approximated 
steady state is anticipated to converge to the exact steady 
state upon increasing the size of the portion of leads in- 
cluded in the projected region. In fact, one can presum- 
ably neglect any direct influence that the interaction in 



the scattering region may give to the leads sufficiently far 
from the impurity because of the screening effect. 

Finally, we observe that, from the technical point of 
view, the complexity of the Gutzwiller problem is the 
same as in the equilibrium case. The only difference is 
that the calculation of the ground state of the quadratic 
Hamiltonian Hq[R, A], see Eq. (|6"Tj) , is replaced with the 
calculation of its steady state. It follows that the techni- 
cal limit of the approach at finite bias is the same as in 
the equilibrium case: the number of variational param- 
eters scales exponentially with the number of orbitals of 
the junction. 

D. Possible connection with the Hershfleld 
effective equilibrium theory 

In his seminal work 5 ^ Hcrshfield proposed an alterna- 
tive description of the steady state in terms of an effective 
equilibrium theory defined by a modified Hamiltonian of 
the form 

H$ =H + . (63) 

It can be proven^ that the nonequilibrium steady state 
can be expressed in a unique way in the Boltzmann form 

p$ oc cxp(— , (64) 

where the so called "bias operator" y encodes the de- 
pendency on the bias. In the latest years the interest in 
this alternative formulation has grown considerably. In 
fact, beyond its purely conceptual relevance in the gen- 
eral context of nonequilibrium statistical mechanics ) 52 ' 53 
several promising numerical techniques have successfully 
implemented this scheme i 30 ' 31 ' 33 ' 34 

In this section we discuss a general aspect of our vari- 
ational scheme in relation with the effective equilibrium 
Hcrshfield framework. 

Let us consider a Gutzwiller wavefunction of the form 

|*a>=PgI*o>, (65) 

where A is a subrcgion of the system that contains the 
impurity, and Vq is the most general operator generated 
by any algebraic combination of local fermionic ladder 
operators c' aRa . and c aRa with fie A. 

As discussed in the previous section, we do not ex- 
pect that the interaction affects directly the electrons 
sufficiently far from the scattering region — because of 
the screening effect. This observation suggests that the 
Gutzwiller steady state variationally determined from the 
ansatz [Eq. |65|) ] 

\* AS )=V G S M S ) (66) 

is a good approximation of the exact steady state \^s) if 
A is sufficiently large. In other words, we expect that the 
Gutzwiller projector Vq S produces important variational 



7 



corrections by increasing A only up to a certain finite 
region A (localized around the correlated impurity), and 
that further increase the variational freedom can provide 
only "marginal" corrections. 

The above conjecture is formulated mathematically by 
the following equations: 



Iim(*As| A\V A s) 

|A|t 
|A|f 



(Vs\A\Vs) VAeA (67) 
VleA (68) 

I[*s] ■ (69) 

The meaning of Eq. (j67|) is that the series of Gutzwillcr 
approximations [Eq. (|66p] "approaches" the exact steady 
state in the limit of infinite A. Eq. (|68p expresses the 
existence of the weak limit of the corresponding series 
of Slater determinants I^q 5 )- Finally, Eq. asserts 
that the exact correlated steady state \^s) belongs to 
the same sector of a proper uncorrelatcd Slater determi- 
nant \^° s }, i.e., that |* s ) and |*|) differ only by local 
modifications, see Sec. IIII Al If proven, this would be an 
interesting observation in relation with the general the- 
ory of infinite systems^ 

From the point of view of the Hershfield theory, the 
verification of our conjecture would suggest that, given 
a general correlated impurity problem, the bias operator 
y can be expressed as 



y~y + sy, 



(70) 



where Sy is a local operator, and y$ is a proper "long- 
range" quadratic operator that depends on the Hamilto- 
nian of the system %. 

Note that to a local difference between two Hamilto- 
nians Ha and T-Lb corresponds a "macroscopic" differ- 
ence between the corresponding effective long-range po- 
tential y$ and y~Q. This is consistent with our observa- 
tion that the sector of steady states of Hamiltonians that 
differ only locally belong to different representations of 
the observables, i.e., they are macroscopically different, 
see SecHIOg 

The physical meaning of Eq. (fT0"| is in line with the 
concept of nonequilibrium Nozieres Fermi liquid theory^ 
in the form proposed in Ref. H3 it should be possible to 
represent the properties of the system in terms of weakly 
interacting quasi-particlcs which, by continuity with the 
non-interacting case, should be regarded as scattering 
states of a properly renormalized quadratic Hamiltonian. 



IV. NUMERICAL RESULTS 

In this section we investigate the reliability of the vari- 
ational ansatz [Eq. employing the single orbital An- 
derson impurity model at half-filling as a prototype. 

We define the hybridization function as 



A(e)=E 



v k 2 



1 



i0+ 



(71) 



and the hybridization width as 



V? 



(72) 



The half-bandwidth W of the leads will be used as the 
unit of energy. 

A. The steady state 



In Sec. IIII 51 we derived the Gutz wilier nonequilibrium 
steady state for the single orbital Anderson impurity 
model. We also demonstrated that in the steady state 
the phase #02 is zero, see Eq. (|4"5)) . and that P02 is deter- 
mined by the condition 



2 VPo2(1-Po 2 ) 







where 



R = 2JP 02 (l-P 02 ). 



(73) 



(74) 



see Eqs. (@5 and JT 
can be shown that 



Using, e.g., scattering theory, it 



(*?|v|*£ 



where 



VZ J a=-l 



Z = R 2 



a 



(75) 



(76) 



can be interpreted as the quasiparticle weight of a single- 
particle cxcitation ) 54 i 55 and 



-Im 



e + i0+ - ZA(e + iO J 



(77) 



has the form of a renormalized uncorrelated impurity 
spectral function. 

It is convenient to summarize here some of the 
Gutzwiller results derived in Ref. H3 for the steady state 
and to underline the limits of the method from the quan- 
titative point of view. This will facilitate the subsequent 
analysis of the dynamics. 

In the rest of this section we will assume the so called 
wide band limit, i.e., V << W. 

At $ = and large U/W one finds that 



W 



cxp 



7T U 

16 r 



Tg 

r 



(78) 



The above expression for Tjf can be interpreted as the 
"Gutzwiller approximation" for the Kondo temperature. 

Notice that Tj^ has two fundamental differences with 
respect to the correct value of the Kondo temperature^ 



/ T f irU 7T r 

T K - U\ — exp 1 

\ 2U y \ 8 r 2 U 



(79) 
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(i) the universal prefactor in the exponent should be 7r/8 
and not 7r/16, and (ii) the factor W/T (which is equal to 
1/r in our units) in Eq. (j78|) diverges in the wide band 
limit. The divergence of for W — > oo reflects the 
unreliability of the method in this limit. In fact, it can 
be verified^ that the Gutzwiller approximation predicts 
that Z — > 1 when W/T — > oo independently on the value 
of U. 

At finite <fr the steady-state quasi-particle weight Z 
vanishes for a finite U — U c . This (unphysical) critical 
point occurs when 



(80) 



that is out of the regime of validity of the method^ 

When the method is extended to out of equilibrium, 
the inaccuracy in Z is reflected in the inaccuracy in the 
steady state current 



.v- Vfc. 



;*i4c fcCT ,_xi*) 



(81) 



In fact, when |^) is a Gutzwiller wavefunction of the 
form [Eq. ([5])] the expectation value of the current in the 
steady state is given by2£ 



deZT(e)pl(e) 



(82) 



On the other hand, the method provides results in ac- 
cordance with the expected universality in terms of the 
Gutzwiller Kondo energy scale T^ defined in Eq. (|78l) . 
In particular, it can be shown^i that the approximated 
Gutzwiller conductance G at zero bias is universal, and 
that its curvature at $ = has the expected Tjf depen- 
dence: 



d 2 G 
d$> 2 



-1 



(83) 



B. Interaction quench 

In this section we study the dynamics of the Anderson 
system [Eq. prepared in its noncquilibrium steady 
state at given bias $ after a sudden change in U at the 
initial time t = 0. This way of perturbing the system is 
commonly referred to as an interaction quench. 

In order to assess the reliability of the variational 
ansatz defined in Eq. ([8]) we focus on interaction quenches 
from U — 0. This is convenient as, for such a non- 
interacting system, the Gutzwiller steady-state formal- 
ism reproduces the exact solution^ In our calculations 
we assume a flat density of states and T/W — 0.1 (essen- 
tially the wide band limit). 

When the interaction is turned on, a disturbance prop- 
agates away from the impurity through the leads with a 
speed that is given by the Fermi velocity of the leads 
VFr- and every local observable converge to its steady 



state value after some relaxation time due to dissipa- 
tion in the leads 4^ We remark that, unlike the case of 
the Hubbard model^S the Gutzwiller parameter P02 - 
which is the expectation value of the local operator P02, 
see Eq. (fTTj). — and its conjugate variable #02 generally 
relax to proper steady state values, see Fig. [TJ 

It can be verified that the curvature of the current 
at t = obtained from Eq. (gS]) and Eqs. d22)-([Ml) is 
quantitatively correct at very short times, as 



Kt 



(*o| J]] |*o) 



(84) 



The fact that the approach provides quantitatively re- 
liable results at short times is not surprising, as the 
Gutzwiller projector is expected to adequately describe 
the physics of Coulomb blockade before the disturbance 
caused by the interaction quench has propagated away 
from the impurity, and other correlation effects has had 
time to build up and affect the dynamics of the system. 

In order to discuss the reliability of our approximation 
also at longer times we have calculated the time evolution 
of the interacting current I(t) divided by its value at 
t = for different biases and interaction strengths, see 
Figs [2] and [3] Our calculations are compared with the 
numerically-exact Monte Carlo results of Ref. [H] in two 
different regimes of parameters: the weak bias regime 
$ < Tk, and the large bias regime, $ > Tk- 

At $ < Tk, sec Fig. [51 the quantum Monte Carlo nor- 
malized current I(t)/I(0) is essentially independent on 
the bias at short times, with a pronounced undershoot 
at tr ~ 1. The Gutzwiller time dependent current is 
in good agreement with the Monte Carlo results up to 
tr < 0.4. At longer times, instead, the Gutzwiller cur- 
rent deviates from the correct solution. In particular, the 
double-occupancy oscillations in the impurity developed 
after the interaction quench are very weakly damped (see 
Scc lIVDj) . and induce unphysical oscillations also in the 
current. Notice that in the weak bias regime the mag- 
netic correlations, that are not properly taken into ac- 
count in our approximation, play a very important role. 
A Kondo cloud builds up and damp the charge fluctua- 
tions in the impurity induced by the interaction quench 
after a time T given roughly by T < 1/Tk- In Fig. [5] 
this corresponds to Tr ~ 2.3 and 4.7 for U/T = 4 and 6 
respectively. Indeed, we interpret the pronounced mini- 
mum of the current at tT ~ 1 and the subsequent small 
oscillations visible at longer times for U/T = 6 in the 
Monte Carlo calculations as a consequence of oscillations 
in the double occupancy. A better variational descrip- 
tion of the magnetic correlations between the impurity 
and the surrounding lead electrons is necessary to repro- 
duce this damping. 

At $ > Tk , see Fig. [3J the the quantum Monte Carlo 
current depends on the bias even at short times. The 
time T required for convergence to the steady state is 
much shorter than in the small bias regime, and is essen- 
tially independent of U~ In this case the Gutzwiller time 
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0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 

Re[P 02 e^] 

FIG. 1: (Color online) Time evolution, for F = 0.1, of the 
variational parameters P02 and #02 in an interaction quench 
from U/T — to U/T = 4.0, starting from the non-interacting 
steady state with bias, $/T = 4.0. P02 starts from the non- 
interacting value Po2(t = 0) = 1/2 and evolves toward the 
interacting steady state value Po2(i — ► 00) « 0.27. 




FIG. 2: (Color online) Interaction quenches from U/T = for 
t < to U/T = 4 and 6 for t > 0. The scaled Gutzwiller tran- 
sient current I(t)/I(0) (solid lines) are shown for two biases 
<E>/r = 1/16 and 3/4 (higher bias giving higher I(t — > 00), in- 
dicated by horizontal lines). For comparison the correspond- 
ing quantum Monte Carlo results from Ref.[I3 are shown (dot- 
ted lines). 



dependent current is in good agreement with the Monte 
Carlo calculations only at very short times, tT < 0.1. For 
longer times the damping rate of the current oscillations 
is underestimated as in the weak bias regime. Note that 
the influence of noncquilibrium effects around the impu- 
rity is strong in this regime. Indeed, it is known that the 
bias can destroy Kondo effect when <!> > Tk ^ We believe 
that the fast damping of the oscillations of the quantum 
Monte Carlo current is a manifestation of this nonequi- 
librium effect, which is sufficiently strong to contrast the 
formation of a magnetic moment already at tT <C 1. Al- 
though the destruction of the Kondo resonance can be 
qualitatively described already in our approximation^ a 
reliable description of the quantum fluctuations induced 
by the bias in the vicinity of the impurity would require 
an extended Gutzwiller projector. 

In conclusion, the Gutzwiller dynamics is quantita- 
tively accurate at sufficiently short times (depending on 
the regime of parameters considered). At longer times, 
instead, the Gutzwiller dynamics is no more reliable. 
In particular, the charge oscillations developed by the 
current after the interaction quench are very weakly 
damped. 



C. Bias quench 

In this section we study the dynamics of the system 
prepared in a given equilibrium configuration after a sud- 
den shift of the chemical potentials fi a in the leads; a 
so called bias quench. We limit the discussion to sym- 
metric quenches from zero to finite chemical potentials 
fi a = a$/2. 

In our calculation we adopt the following "smoothed" 




0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 

tT 



FIG. 3: (Color online) Interaction quench from U = to 
U/T = 4 and 6 (counting from above) at a bias <J>/r = 4 and 
r = 0.1 (solid lines). The corresponding steady-state currents 
are indicated by horizontal lines. For comparison, quantum 
Monte Carlo results from Ref . Qi] for the same system (dotted 
lines) . 

square density of states in the leads 

= 2 1 + e^-^' 1 + e - v ^ +w "> 

with vT = 3, where v is a measure of smoothness of the 
edges of the DOS and W = 1 is the half-bandwidth. 

In our time dependent calculations the system is ini- 
tially prepared in its approximated interacting Gutzwiller 
equilibrium state. Unfortunately, unlike the case of the 
interaction quench, the initial state can not be prepared 
exactly in the Gutzwiller approximation, as U is finite 
even at t = 0. For this reason we can not attempt a 
quantitative verification of the Gutzwiller dynamics in 
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FIG. 4: (Color online) Current / for bias quenches from the 
unbiased state $/r = to finite bias <£>/r = 0.25, 0.5 and 
1.0 (higher bias correspond to higher current) for U/T = 4. 
The Gutzwiller results (solid lines) are compared to the non- 
interacting quench (dash dotted lines) and quantum Monte 
Carlo results from Ref. [T^ (dotted lines) . 
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FIG. 5: (Color online) Steady state currents I s in the weak 
bias regime. The Gutzwiller results (circles) are compared 
to the non-interacting bias quench (squares) and quantum 
Monte Carlo results (diamonds) from Ref. [l(J. The effect of 
using the "a priori" estimate of the Kondo temperature to fix 
R is also shown (triangles). 



this case, not even at short times. The transient current 
after the application of weak biases $ < T is shown in 
Fig.|l]for three different interaction strengths, U/T = 0, 4 
and 6. In this weak bias regime the Gutzwiller parame- 
ters P02 and #02 do not change appreciably after the bias 
quench, independently of T and U. Consequently, the 
current evolves in time essentially as in a renormalizcd 
uncorrelated system W§, see Eq. (fl~5]) . with a constant 
renormalization R giving an overall reduction in the cur- 
rent at finite U . In the special case of U = the renor- 
malization constant R is equal to 1, and the Gutzwiller 
calculation reproduces the exact dynamics of the non- 
interacting system. 

The shoulder like features in the transient current in 
Fig. |4j that are present even in absence of interaction, 
are due to resonance effects induced by the specific form 



of the density of states^ and would not be visible in the 
wide band limit . 

From the comparison of our calculations with the quan- 
tum Monte Carlo results of Ref. [l^ it is clear that the 
suppression of the current due to U is severely underes- 
timated in the Gutzwiller approximation. This is due to 
the problem discussed in Sec. IIV Al that R ~ 1 for large 
bandwidths, even when U > T, causing the Gutzwiller 
dynamics to be essentially generated by Hq for weak bias 
quenches. 

Although a quantitative variational description of the 
dynamics is impossible with the oversimplified ansatz 
[Eq. (JSJ], the fact that the Gutzwiller current behaves 
essentially as for a renormalized uncorrelated system H§ 
for small bias quenches suggests a simple interpretation 
of some of the qualitative features of the current provided 
by the Monte Carlo calculations. In Fig. [3] is shown the 
steady state current generated after a bias quench by T~L^, 
with 



R = \/Tk/T , 



(86) 



where Tk is given by Eq. (fT^l) rather than the Gutzwiller 
approximation. This simple calculation better repro- 
duces the steady state current of the quantum Monte 
Carlo calculations. This supports our statement that the 
main problem of the time dependent Gutzwiller method 
for small bias quenches is that the renormalization factor 
R is underestimated in the initial equilibrium state. 



D. Note about relaxation and long-time behavior 

As we have shown, the Gutzwiller dynamics is able to 
capture dissipation in the leads, and a Gutzwiller steady- 
state can generally be reached after a sufficient thermal- 
ization time. In this section we discuss the long-time 
dynamics more in detail, and we analyze its dependency 
on the tight binding parameters of the Hamiltonian and 
the bias. 

In Sec. IIV Al we anticipated that at finite bias the 
steady-state theory is qualitatively reliable only for U <C 
U c , where U c is the spurious critical interaction discussed 
in Sec. IIV Al characterized by a vanishing quasi-particle 
weight Z = R 2 . For this reason we confine our analysis 
to U < t/ c , which corresponds to the weak bias regime 
defined by the condition 



(87) 



see Eq. ([50]). 

Within the range of parameters considered, we find 
that the long-time behavior of Po2(t) is, to a good ap- 
proximation, an exponentially damped oscillation with a 
measurable frequency / and a relaxation time r. 

Interestingly, we find that both r and / are indepen- 
dent on the initial condition, and depend solely on the 
bias and the tight binding parameters of the Hamilto- 
nian (not shown). The dependency of r and / on U is 
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FIG. 6: (Color online) Upper panel. Oscillation frequency 
/ and inverse relaxation time 1/r as a function of UT for 
zero and finite bias, $/r — and 1 (full and shaded lines 
respectively). The dotted line represents the frequency fo 
defined in Eq. (|91|l . Lower panel. Steady state values of R 
and V. The maximum of r _1 at $ = 0, that is realized at 
UT ss 12.8, is indicated by the vertical line. 



illustrated in the upper panel of Fig. [5] for zero and finite 
bias. 

Our numerical results show that the Gutzwiller dy- 
namics spontaneously reaches a steady state except in 
the following three limiting cases: (i) at U = 0, (ii) at 
U —> oo and $ = 0, and (hi) for U > U c at finite bias. 
The inverse relaxation time r _1 increases by increasing 
U from up to a maximum value, and vanishes after- 
wards. Note that the effect of the bias is to decrease the 
relaxation time for any given value of U. The frequency 
of the oscillations / grows monotonically with U, and is 
smaller at finite bias with respect to the equilibrium case. 

In order to better discuss our numerical results we lin- 
earize Eqs. (f2"Tj) - ([2"4"| around the steady state solution^ 
We introduce steady state values (indicated by bars) and 
deviations from the steady state 

P02 = P02 + AP02 

#02 = #02 + A6* 2 



V = V + AV 
It can be easily shown that 

d t AP 02 = RVA6 02 
d t A6 2 = 



(88) 



U 4V 
2V R 3 



(89) 



d?AP Q[ 



UR 



-AV 



4V^ 



■AP„ 



(90) 



Notice that if UR = Eq. J90|) reduces to the differ- 
ential equation of a simple harmonic oscillator for AP02 , 
with a resonance frequency 



, _ 1 2|V| 
/o= 2^~fT 



(91) 



This observation explains why the expectation value of 
the double occupancy does not relax whenever either U 
or R is zero. In particular, Eq. ([9"0"]) shows that the crit- 
ical point at finite bias, see Eq. (|78|) . is reflected by a 
spurious dynamical transition, characterized by a diverg- 
ing relaxation time. 

The result that the relaxation time of the system di- 
verges in the large U limit at $ = is physically cor- 
rect, as in this limit the only energy scale is the Kondo 
temperature, which vanishes exponentially with U, see 
Eq. (|79[) . Also the fact that the relaxation time r de- 
creases by increasing the bias^, which reflects the fact 
that bias "destroys" Kondo effect r& is qualitatively cap- 
tured by the Gutzwiller approximation. On the contrary, 
the fact that the expectation value of the double occu- 
pancy is unable to relax at U = is clearly a drawback 
of our oversimplified variational ansatz4^ This incorrect 
feature of the Gutzwiller solution questions the whole 
observed behavior of t" 1 for UT < 12.8 in Fig. [5] 

Let us now consider the behavior of the frequency /. 
From Fig. [S]we see that / ~ /o for all values of U. Using 
Eq. (|75|) . it can be readily shown that 



hm ( ±m 

U^oc V 2-7T R 



00 . 



(92) 



Eliminating A#o2 we obtain the following equation for 



This limit exposes another drawback of the approxima- 
tion, as the same dimensional argument used above for 
the relaxation time should be applicable to the period of 
oscillation: at large U the only time scale is the inverse 
of the Kondo temperature, which vanishes exponentially 
for U — > 00. For this reason we conclude that the be- 
havior of / can not be addressed within our variational 
approximation. 

In conclusion, the simple variational ansatz [Eq. ([5])] 
is able to capture dissipation in the leads, but a more 
general type of variational wavefunction is necessary in 
order to obtain a reliable description of the long time 
behavior. 



V. CONCLUSIONS 

We have generalized the time dependent Gutzwiller ap- 
proach of Ref. HH to impurity problems, modeling a quan- 
tum system, e.g., a molecule or a quantum dot, coupled 
to metallic leads. The dissipation effects in the junction 
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due to the coupling to the leads are naturally accounted 
for by the method, and a steady state is reached spon- 
taneously, without, e.g., including any fictitious bosonic 
bath^ 

We have shown that the time dependent theory, that 
is based 4 ^ on the Dirac-Frenkel variational principle ; 4 ^ 41 
enables a natural derivation of a variational principle for 
the steady state, that reproduces the ansatz previously 
proposed on different grounds in Ref. H3- We believe that 
our formulation of the steady state problem represents 
an important conceptual advancement in itself, and that 
the basic idea behind it could be exploited also to develop 
different numerical techniques. 

The Gutzwiller method for the steady state is partic- 
ularly appealing from the computational point of view 
as, although there are important conceptual differences 
between the Gutzwiller theory for the nonequilibrium 
steady state and the conventional method at equilibrium, 
it seems to us that these differences are not accompanied 
by any substantial increase in computational complexity 
— which is generally the case for other methods. 

We have investigated the reliability of the time depen- 
dent variational method in the case of a Gutzwiller pro- 
jector acting only on the impurity, see Eq. (|5]l. which is 
the simplest conceivable type of Gutzwiller variational 
function. From the comparison with the quantum Monte 
Carlo data of Ref. [l^ we concluded that the obtained 
Gutzwiller dynamics after an interaction quench is quan- 
titatively accurate at sufficiently short times. Neverthe- 
less, due to the limited description of the quantum fluc- 
tuations and magnetic correlations between the impurity 



and the surrounding electrons, the damping rate of the 
current oscillations generated in the impurity by the sud- 
den rise of the interaction U at t = is underestimated 
in this approximation, and the steady state current is 
quantitatively incorrect. 

Considering the fact that the disturbance due to the in- 
teraction quench propagates with a finite velocity v ~ vf 
from the impurity, the behavior of the current will be bet- 
ter described, at least at short times, by extending the 
region of action of the Gutzwiller projector to a portion 
of the leads surrounding the impurity. But we also ex- 
pect that this more elaborate variational ansatz would 
provide a better description of the time dependent cur- 
rent at any time, uniformly. In fact, the correction due 
to the Gutzwiller projector on the correlation functions 
is expected to vanish at distances sufficiently large from 
the impurity due to the screening of the Coulomb interac- 
tion provided by the conduction electrons in the vicinity 
of the impurity. The numerical verification of this hy- 
pothesis will constitute an interesting extension of this 
work. 
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